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Melanophryniscus is a bufonid frog genus with a broad geographic distribution over southeastern South America. 
In recent years, several new species of Melanophryniscus have been discovered in southern Brazil showing a 
distinctive life-history strategy for the genus - breeding in phytotelmata - as well as a strong association with 
high-altitude regions. In this study, we use mitogenomic data to infer the phylogenetic relationships among 
diurnal, phytotelm-breeding Melanophryniscus and to determine the timing of their divergence. We obtained the 
mitochondrial genomes (not including the control region) for eight individuals of Melanophryniscus representing 
all three described species (M. alipioi, M. milanoi, and M. xanthostomus), as well as some recently-discovered and 
potentially new species. Gene order was conserved in all species and corresponded to the general order found in 
bufonids. Although the phylogenetic relationships among the studied species was poorly supported, dating 
confirmed that they diverged during the Pleistocene, suggesting that phytotelm breeding could have arisen 
during drier periods in the glacial/interglacial cycles due to a decrease in the availability of permanent streams 


or ephemeral/temporary streams or ponds in which Melanophryniscus species commonly breed. 


1. Introduction 


Melanophryniscus (Anura: Bufonidae) is a fascinating frog genus 
with a broad geographic distribution over southeastern South America. 
With 29 currently described species, the genus is distributed across the 
south and southeastern Brazilian Atlantic Forest, wetland and grassland 
regions of Brazil to the inter-Andean valleys in Bolivia, and areas across 
Paraguay and Uruguay down to central Argentina (Frost, 2016). The 
genus is a charismatic component of the anurofauna of the Neotropics, 
with many species with bright colours that advertise the alkaloid-con- 
taining defensive chemicals secreted by the granular glands in their 
skin, which are either sequestered from the arthropod in their diet or 
endogenously biosynthesized (Hantak et al., 2013; Jeckel et al., 2015). 
Although many species of Melanophryniscus are still poorly studied, 
there has been increasing concern regarding the conservation of its 
species. For instance, of all 23 species of Melanophryniscus currently 
assessed on the IUCN database, 11 are under some level of threat 


(IUCN, 2016). 

Although most Melanophryniscus species breed either in small per- 
manent streams or on ephemeral/temporary streams or ponds (see 
Baldo et al., 2014), it was recently discovered that some species evolved 
in association with plants for phytotelm breeding (Langone et al., 2008; 
Steinbach-Padilha, 2008; Bornschein et al., 2015). As a consequence, 
the distribution of phytotelm-breeding species is not constrained by the 
presence of “larger” bodies of water, thus allowing these species to 
occupy other types of habitat, such as campos de altitude (highland 
grasslands), in which there are no sources of accumulated water other 
than phytotelms (Langone et al., 2008). There are five currently re- 
cognized species of phytotelm-breeding Melanophryniscus, all which 
from the highlands of the Atlantic Forest of the states of Parana and 
Santa Catarina, southern Brazil, namely M. alipioi (Langone et al., 
2008), M. vilavelhensis (Steinbach-Padilha, 2008), and M. biancae, M. 
milanoi, and M. xanthostomus (Bornschein et al., 2015). Based on pre- 
liminary data on the phylogeny of the genus (Baldo et al., 2014), all five 
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phytotelm-breeding species and the pond-breeding M. moreirae form a 
monophyletic group. Moreover, Firkowski et al. (2016) showed that M. 
biancae is phylogenetically close to M. vilavelhensis and that this group is 
phylogenetically distant from M. alipioi, M. milanoi, and M. xanthos- 
tomus. Although the phylogeny of Firkowski et al. (2016) does not in- 
clude M. moreirae, we tentatively propose that M. biancae and M. vila- 
velhensis should be considered sister species given that, in addition to 
their phylogenetic proximity, both species live in the same habitat 
(grasslands; see below), breed in Eriocaulon ligulatum, seem to be noc- 
turnal, and have similar morphology (e.g. small snout-vent length and 
general colour pattern) (Steinbach-Padilha, 2008; Bornschein et al., 
2015; MRB per. obs. regarding the species used M. vilavelhensis during 
the breeding season). On the other hand, the remaining phytotelm- 
breeding Melanophryniscus are diurnal, have larger body size, and re- 
produce in bromeliads (on rare occasions also on fallen dead bamboo), 
both in forests and in grasslands (Langone et al., 2008; Bornschein 
et al., 2015). It is important to note that kind of habitat in the grass- 
lands with marshes where that pair of species occurs is distinct from 
that in the grasslands at the top of the mountains that harbour M. alipioi 
(Langone et al., 2008) and some other newly discovered related po- 
pulations, which are classified as “Refúgio Vegetacional”, whereas M. 
biancae and M. vilavelhensis are found in “Estepe Gramineo Lenhosa” 
(sensu Veloso et al., 1991). Although the fact that these different clas- 
sifications do not necessarily imply that they have a different phyto- 
geographic origin, the “Refigios Vegetacionais” of the mountains of the 
Serra do Mar of Parana and Santa Catarina do not contain the marshes 
with E. ligulatum where M. biancae and M. vilavelhensis reproduce (MRB 
pers. obs.). Therefore, based on ecological and morphological traits, one 
can define two groups phytotelm-breeding of species of Melano- 
phryniscus: the nocturnal M. biancae and M. vilavelhensis, which breed 
on Eriocaulon, and the diurnal M. alipioi, M. milanoi, and M. xanthos- 
tomus, which breed on bromeliads. However, based on the results of 
Firkowski et al. (2016), it is possible that these three species are in fact 
complexes of highly endemic species. 

An important step towards understanding diversification in 
Melanophryniscus is to assess the timing of recent speciation events, so 
that they can be interpreted in relation to alternative divergence me- 
chanisms. In this study, we investigate the evolution of diurnal, phy- 
totelm-breeding species of Melanophryniscus using mitogenomic data. In 
particular, our goals are threefold: (1) to describe the structure of the 
mitochondrial genome of phytotelm-breeding Melanophryniscus; (2) to 
infer the phylogenetic relationships among populations of the three 
recently described species of Melanophryniscus, as well as some poten- 
tially undescribed new species, based on their mitogenomes; and (3) to 
infer the timing of their diversification. 


2. Methods 


Tissue samples were obtained from field-collected specimens in 
seven locations for Melanophryniscus, including both described and 
potentially undescribed species (Table 1). Given that, although our 
datasets involve mitogenomes, they only involve a single locus and 
therefore is limited regarding species delimitations. However, given the 
evidence in Firkowski et al. (2016), we tentatively called them either by 
the available name when the sample was obtained from the type lo- 
cality, or as “affinis”, when data from Firkowski et al. (2016) indicated 
that it could be a distinct species. Voucher specimens were deposited in 
the herpetological collection of the Department of Zoology of the Uni- 
versidade Federal do Paraná (DZUP) and Museu de História Natural 
Capão da Imbuia (MHNCI), Curitiba, Brazil. Total genomic DNA was 
extracted using PureLink™ Genomic DNA kit (Invitrogen™, USA), ac- 
cording to the manufacturer's instructions. mtDNA sequences were 
obtained as off-target regions from another study (Pie et al., un- 
published results) using target capture of Ultraconserved Elements (see 
Faircloth et al., 2012). We used the resulting reads from the assemble 
performed by Trinity (Grabherr et al., 2011) and executed in Phyluce 
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(Faircloth, 2013). As expected (Hung et al., 2013), the longest contig 
included the entire mtDNA genome - except for the control region. 
mtDNA genome notation was carried out using Miros (Bernt et al., 
2013). Sequence reads for this project are available from NCBI Bio- 
Project PRJNA391191. 

We obtained additional mtDNA genome sequences from GenBank 
for the only other Melanophryniscus mitogenome available to date (M. 
simplex, Machado et al., 2016), and six other bufonid species (four Bufo 
species, one Duttaphrynus, and one Leptophryne), as well as 10 species of 
Nobleobatrachia (as defined by Frost et al. (2006:196)), representing 
other six families (Table 1). Sequence alignment of coding genes was 
carried out using Marrt v. 5 (Katoh et al., 2005) using the TRANsLaToRX 
server (Abascal et al., 2010), which first translated sequences into 
aminoacids prior to the alignment and back-translated into nucleotides 
for later analyses. Non-coding genes (rrnS, rrnL, tRNAs) were aligned 
directly in Marrr 7.2 (Katoh and Standley, 2013). The best partitioning 
scheme was then determined using PARTITIONFINDER v. 2 (Lanfear et al., 
2016) based on its greedy algorithm (Lanfear et al., 2012), with models 
selected according to the Bayesian Information Criterion (Schwarz, 
1978). We performed different analyses to choose the models to be used 
in tree estimation methods by Bayesian inference (BI) (constrained to 
the 24 models of evolution available in MrBayes or models available in 
Beast v. 1.8) and maximum likelihood (ML) (constrained to the GTR + T 
model for RAxML). 

We estimated the phylogeny using RAxML v. 8.2.4 (Stamatakis, 
2014), using the selected partitioning scheme for the ML analysis, and 
accessed nodal support with 1000 bootstrap replicates. We inferred the 
phylogenetic relationships among the studied species using BI in 
MrBayes v. 3.2.6 (Ronquist et al., 2012), selecting the partitioning 
scheme and models selected with PartiTionFinpER (with two independent 
replicates of 10” generations, with one cold and seven heated chains, 
with sampling every 1000th generations). We examined the con- 
vergence of both runs in Tracer v. 1.6 (Rambaut and Drummond, 2013; 
Rambaut et al., 2013) and Awry (Wilgenbusch et al., 2004; Nylander 
et al., 2008). As both runs appeared to converge, we combined the trees 
of both replicates after discarding the first 25% samples of each, as 
burn-in. 

To test if the same molecular clock model could be implemented for 
all the partitions defined for the Bayesian analyses, we estimated ML 
trees in RAxML under the GTR +T model for all the defined subsets, 
and provided these trees as input to estimate a minimum branch-score 
distance matrix and the optimal number of clock partitions using 
CiocksTaR2 (Duchéne et al., 2014) in R v. 3.3 (R Development Core 
Team, 2016). The timing of lineage divergence was then inferred using 
Beast v. 1.8 (Drummond and Rambaut, 2007; Drummond et al., 2012), 
with two independent runs of 50 x 10° generations with a Yule prior, 
sampled every 1000th generation, with the first 25% omitted as burn- 
in. Given that there are no known Melanophryniscus fossils, we em- 
ployed a relaxed log-normal molecular clock (Drummond et al., 2006) 
to calibrate the age of the split between Melanophryniscus and the re- 
maining Bufonidae, set as a normal prior with a mean age of 
67.92 million years ago (Mya) and distribution between 52.7 and 
92.7 Mya (see Van Bocxlaer et al., 2010; Firkowski et al., 2016). Con- 
vergence of the runs was examined in Tracer v. 1.6. As the runs with the 
complete dataset failed to converge, possible due unstable topology at 
the root of the tree, we performed the dating analysis including only the 
Bufonidae. As the two runs with the Bufonidae-only dataset converged, 
we combined both runs (after removal of a 25% burn-in), and analyzed 
the combined dataset to obtain a majority rule consensus tree, using 
TREEANNOTATOR V. 1.8 (Drummond et al., 2012). We then used DensitTREE 
v. 2.2 (Bouckaert, 2010; Bouckaert and Heled, 2014) to explore the 
differences between the topologies obtained in the different analyses. 
We performed all phylogenetic analyses at the Cres server (Miller 
et al., 2010). 
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Sampled taxa used for phylogenetic analyses and corresponding GenBank accession numbers. Localities are provided for the species sequenced in the present study. 


Taxa Locality Family GenBank accession 
Melanophryniscus sp. aff. alipioi DZUP345 Apiai, Sao Paulo Bufonidae KY260875 
Melanophryniscus sp. aff. alipioi DZUP321 Itapiroca, Parana Bufonidae KY260872 
Melanophryniscus alipioi DZUP344 Capivari Grande, Parana Bufonidae KY260870 
Melanophryniscus sp. aff. alipioi DZUP201 Torre da Prata, Parana Bufonidae KY260871 
Melanophryniscus sp. aff. alipioi DZUP323 Morro dos Padres, Parana Bufonidae KY260873 
Melanophryniscus sp. aff. xanthostomus MHNCI 9806 Condominio Vale dos Lagos, Santa Catarina Bufonidae KY260876 
Melanophryniscus milanoi DZUP200 Morro do Bati, Santa Catarina Bufonidae KY260874 
Melanophryniscus sp. aff. milanoi DZUP437 Morro do Cachorro, Santa Catarina Bufonidae KY260877 
Melanophryniscus simplex - Bufonidae KT221611 
Bufo gargarizans - Bufonidae NC008410 
Bufo japonicus - Bufonidae NC009886 
Bufo stejnegeri - Bufonidae KR136211 
Bufo tibetanus - Bufonidae NC020048 
Duttaphrynus melanostictus - Bufonidae NC005794 
Leptophrynine borbonica - Bufonidae JX564876 
Espadarana prosoblepon - Centrolenidae JX564857 
Hyalinobatrachium fleischmanni - Centrolenidae JX564869 
Ceratophrys ornata - Ceratophrydae JX564858 
Telmatobius bolivianus - Telmatobidae NC020002 
Rhinoderma darwinii - Rhinodermatidae JX564891 
Hyla japonica - Hylidae NC010232 
Leptodactylus melanonotus - Leptodactylidae JX564873 
Phyllomedusa tomopterna - Hylidae JX564887 
Dendrobates auratus - Dendrobatidae JX564862 
Mannophryne trinitatis - Dendrobatidae JX564878 


3. Results 


We successfully obtained eight mitogenomes, varying between 
15,206 and 15,210 base pairs without the control region (Table 1), with 
an average coverage of 30.7 (range = 8.44-65.82). These genome sizes 
are shorter than in the African clawed frog Xenopus laevis (c. 17,500 bp, 
GenBank accession NC001573) due to the control region not being 
computed. All the typical 37 genes are present (Fig. 1) and the overall 
nucleotide composition is biased against guanine and cytosine 
(G + C = 38.9%). Only the nad6 gene and eight tRNAs (tmE, trnS2, 
trnY, trnC, trnN, trnA, trnQ and trnP) are encoded in the heavy strand. 
The trnV gene separates the ribosomal RNAs (rrnS and rmL). The length 
of rrnS is the same for the three species, 934 bp, while the length of rrnL 
slightly fluctuate from 1596 to 1599 bp. Both RNAs sequences are A 


trnR 


Fig. 1. The structure of the mitochondrial genome of diurnal, phytotelm-breeding 
Melanophryniscus. 


(adenine) and T (thymine) rich with proportions around 56.1% for rrnS 
and 63.5% for rmL. 

The aligned dataset had 15,319 bp. The partitioning analyses rea- 
lized in PartitionFinder returned a scheme with 13 partitions for 
MrBaxes (Table S1), 14 for Beast (Table S2) (10 for the Bufonidae-only 
dataset; Table S3), and 11 for RAxML (Table S4). The phylogenies in- 
ferred by BI and ML using the nucleotide alignment presented similar 
results (Figs. 2 and 3), with Melanophryniscus being recovered as sister 
to other analyzed bufonids (Bufo, Duttaphrynus, and Leptophryne). Both 
analyses recovered the family Centrolenidae (represented by Hyalino- 
batrachium and Espadarana) as sister to the Bufonidae, but this re- 
lationship was not well supported in the ML tree. A sister group re- 
lationship between Telmatobidae and Ceratophrydae, and their close 
relationship with Rhinodermatidae were recovered by both methods 
with high support. The two analyses recovered Dendrobatidae as sister 
to a group including Hylidae and Leptodactylidae, but this topology 
was not well supported in the ML phylogeny. Both methods recovered 
M. simplex as sister to all other sampled Melanophryniscus. All phylo- 
genetic analyses place M. alipioi and M. spp. aff. alipioi as most closely 
related to M. milanoi and M. sp. aff. milanoi, whereas M. sp. aff. xan- 
thostomus would be their sister lineage (Figs. 2 and 3). 

The CrocksraR2 analysis returned an optimal number of clock par- 
titions of one, as would be expected for mitochondrial sequences (Fig. 
S1). The dated tree we obtained using Beast for the Bufonidae-only 
dataset showed a different topology than that obtained by the analysis 
in MrBayes and RAxML (Fig. 4A-B). The topology obtained by MrBayes 
and RAxML represented the second most common topology in the Beast 
analysis (Fig. 4C-D). The dated tree returned a Late Cretaceous age for 
the Bufonidae (median 67.99 Mya, 95% HPD 52.74-85.06 Mya; 
Fig. 4A), with a Miocene split between M. simplex and the sampled 
phytotelm Melanophryniscus (median 21.5Mya, 95% HPD 
33.6-12.75 Mya). A much younger age was recovered for the diversi- 
fication of the phytotelm Melanophryniscus (Fig. 4C), estimated to have 
started in the Pleistocene (median 1.68 Mya, 95% HPD 0.95-2.76 Mya), 
whereas the split between M. alipioi + M. sp. aff. alipioi and M. sp. aff. 
xanthostomus was also dated as of Pleistocene age (median 1.31 Mya, 
95% HPD 0.74-2.16 Mya). Our results suggest that the diversification 
of Bufo was already ongoing during the Late Miocene (median 
9.02 Mya, 95% HPD 5.19-14.39 Mya). 
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Fig. 2. Phylogenetic relationships among the studied mitogenomes based on the consensus of two Bayesian inference runs in MrBayes v.3.2.4. of 10” generations sampled each 1000th 


(burn-in of 25%). 
4. Discussion 


The genome content and gene order are consistent among the ex- 
amined Melanophryniscus species, and with other Bufonidae and 
Nobleobatrachia lineages (see Xia et al., 2014). Several studies have 
recovered Dendrobatidae as the sister group of Bufonidae, in general 
with low support (e.g. Hedges and Maxson, 1993; Frost et al., 2006; 
Bossuyt and Roelants, 2009; Pyron and Wiens, 2011). Other authors 
inferred a sister relationship between Bufonidae and Hylidae (Irisarri 
et al., 2010; Kakehashi et al., 2013; Li et al., 2014; Xia et al., 2014; Carr 
et al., 2015). Our results using the nucleotide mitochondrial dataset 
suggests that the sister group of Bufonidae is the family Centrolenidae, 
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as also recovered by Zhang et al. (2013). The sister group relationship 
between Ceratophryidae and Telmatobidae recovered in the present 
study is consistent with the founding of previously published studies 
(e.g. Frost et al., 2006; Bossuyt and Roelants, 2009; Zhang et al., 2013). 
However, it is important to emphasize that inferences regarding the 
relationships among the studied families should be done with caution, 
given the low level of taxon sampling in our analyses. 

As expected, given the calibration used in the present study, the 
estimated Paleocene age for the split between Melanophryniscus and the 
other Bufonidae is similar to those obtained in studies that used similar 
calibrations (e.g. Maciel et al., 2010; Van Bocxlaer et al., 2010), but 
younger than the Late Cretaceous ages obtained by other authors that 
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Fig. 3. Phylogenetic relationships among the studied mitogenomes based on the majority rule consensus of 1000 maximum likelihood bootstrap replicates using RAXML v.8.2.6. 
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Melanophryniscus sp. aff. alipioi DZUP321 


Melanophryniscus alipioi DZUP344 


Melanophryniscus sp. aff. alipioi DZUP345 
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Fig. 4. A. Dated BEAST chronograms of the studied mitogenomes. Node values are of posterior probability; blue bars indicate the 95% HPD for the node age estimates. B. Representation 
of the set of 15,001 postburn-in trees evidencing the three most common topologies (first most common - green, second most common - orange, third most common - black). C. Detail of 
the topology for the Melanophryniscus species from “A”. D. Detail of the set of 15,001 postburn-in trees, from “B”, evidencing the three most common topologies for the sampled 
Melanophryniscus species. Colours are the same as in “B”. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 


used different calibration points (Pramuk et al., 2008; Fouquet et al., 
2012; Portik and Papenfuss, 2015). Such deep divergence between 
Melanophryniscus and other bufonids could be compared to those found 
among some other anuran families (see Bossuyt and Roelants, 2009; 
Zhang et al., 2013), with estimates of an Oligocene age (c. 33 Mya) for 
the onset of the Melanophryniscus diversification (Pyron, 2014), with 
such old divergence and diversification suggesting the possibility that 
this clade could represent a Family level taxon. Our results also re- 
covered a Pleistocene age for the splits between the sampled popula- 
tions of phytotelm Melanophryniscus, in agreement with the results of 
Firkowski et al. (2016), despite the different datasets and methods used 
in our study (concatenated analysis of mitogenomes) and that of 
Firkowski et al. (2016) (three mitochondrial and three nuclear genes, 
species tree analysis), with both using the same calibration. 

Within Melanophryniscus, our results were partially consistent with 
the species tree topology recovered by Firkowski et al. (2016). Our 
MrBayes and RAxML analyses recovered the topology in which M. 
xanthostomus is the earliest-diverging species ((M. alipioi, M. milanoi), 
M. xanthostomus), but with low support, as in Firkowski et al. (2016). 
On the other hand, our Bzasr analysis recovered a topology where M. 
milanoi, the southernmost diurnal, phytotelm-breeding species, was the 
first to diverge, but also with low support (the first mentioned topology 
representing the second most common tree in the Berasr analysis). Al- 
though it is possible that these topological differences were due to the 
different partitioning schemes used for each inference method, the 
analyses in MrBayes and RAxML used different partitioning schemes and 
returned the same topology. Therefore, it seems probable that the 
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differences in the obtained topologies in our study could be related to 
tree rooting, most specifically to the exclusion of the non-bufonid taxa 
for the dating analysis. 

Such young divergence between the sampled populations, within a 
small geographical scale (see map in Bornschein et al., 2015) seems 
consistent with models that propose that endemism patterns in 
southern part of the Atlantic Forest are related to climatic heterogeneity 
(Carnaval et al., 2014). It is possible that Melanophryniscus evolved 
phytotelm-breeding during the drier periods in the Pleistocene cycles of 
glacial and interglacial climate (see Bowen, 2009 and Langone et al., 
2008), a period in which one would expect a decrease in the availability 
of permanent streams or ephemeral/temporary streams or ponds where 
the genus commonly breeds. Likewise, Langone et al. (2008) suggested 
that the reduction in the area of suitable habitat due to the past climate 
change could have isolated populations and promoted speciation, but 
also believe that the isolation of populations in patch of habitat with 
inadequate topography for water to accumulate on the ground forced 
the emergence of reproduction in phytotelmata. 

Supplementary data to this article can be found online at http://dx. 
doi.org/10.1016/j.gene.2017.07.048. 
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